Supporting Information 

Chad M. Topaz 1 '*, Maria R. D'Orsogna 2 , Leah Edelstein-Keshet 3 Andrew J. Bcrnoff 4 

1 Department of Mathematics, Statistics, and Computer Science, Macalester College, 
Saint Paul, Minnesota, United States of America 

2 Department of Mathematics, California State University at Northridge, Los Angeles, 
California, United States of America 

3 Department of Mathematics, University of British Columbia, Vancouver, British 
Columbia, Canada 

4 Department of Mathematics, Harvey Mudd College, Claremont, California, United 
States of America 

* E-mail: ctopaz@macalester.edu 



Model equations 

For convenience, we reintroduce our model equations. Consider a two-dimensional domain fl with spatial 
coordinate x = (x, y). Define p(x, t) = s(x, t) + g(x, t) as the locust population density field, with s(x, t) 
and g(x,i) the solitary and gregarious components, respectively. The locust populations move with 
velocities v Si9 (x, i) and obey the equations 

s + V- (v,*)= -h{p)a + h{p)g, v s = -V(Q s *p), (la) 
g + V • (v g g)= f2(j>)s-fi(p)g, v 9 = -V(Q ff * P ), (lb) 

These equations generalize the classic swarming model 

p t + V • (pv) = 0, v= - / VQ(x-x')p(x',i)dx', (2) 

Jn 

which describes a single population density field advected by a velocity field arising from social interac- 
tions. Eq. ^ has been studied extensively in one and two spatial dimensions for various social interaction 
functions represented by Q, whose negative gradient is the effective social force [?,?,?,?]. Depending on 
Q, solutions include steady swarms, spreading populations, and contracting groups (i.e., blow-up) [?,?,?]. 
In our two-phase model Eqs. {TJ), the velocities are 

v,, ff (x, t) = -VQ s . g *p=- I VQ s , 9 (x - x>(x', t) rfx', (3) 

Jn 

and the social interaction potentials Q s ,g are 

Q s (x-x') = i? s c-l x - x 'l/ r % Qg(x-x') = i? g c-l x - x 'l/' r " -^ g c-l x - x 'l/ Q '. (4) 

Here, R s , R g , A g are interaction magnitudes and r s ,r g and a g are interaction length scales. We require 
RgCLg — A g r g > and A g a 2 g — R g r^ > so that Q g includes short range repulsion and long range attraction, 
as in [?,?,?], as this is the clumping regime, appropriate to capture the tendency of gregarious locusts to 
aggregate. We model the density-dependent rates of interconversion of the solitary and gregarious forms 
as 

1 + [p/ki) 1 + (p/k 2 ) 

The parameters 5i t 2 are maximal rates and k\ t 2 are characteristic locust densities at which the transitions 
occur at half of their maximal values. To the best of our knowledge, our work is the first to consider 
locust phase changes via continuum modeling of locust density [?,?,?,?]. 
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Parameter selection and estimation 

As discussed in the main text, for our numerical results, we use two different sets of phase change 
parameters. For both sets, we use the same social interactions parameters, and we now describe our 
choices for these. 

To estimate R s , R g , and A g , we use explicit velocity computations. The speed of a locust when it is 
alone varies between 72 - 216 m/hr, while the speed of a locust in a group varies in a tighter range of 
144 - 216 m/hr [?]. To make a rough estimate of R 8 , we imagine a hypothetical semi- infinite density field 
p{x, y) = p g rouptl(x) where H(x) is the Heaviside function and, as mentioned in the main text, p gr0 up = 
65 locusts /m 2 is the approximate critical density of a gregarious group [?]. A solitary locust placed at 
the swarm's edge (at the origin) should move to the left with maximal velocity v™ ax = —216 m/hr. From 
Eqn.©, 

v s (Q,0) = {-VQ S * Pgroup n{x)} \ {m = < iax , (6) 

which we solve to find R s = 11.87 m 3 /(hr ■ locust). Similarly, a gregarious locust at the origin should 
move to the right with maximal velocity y™ ax = 216 m/hr, so 

v g (0,0) = {-VQ 3 * Pgroup U(x)} \ m = vf<**. (7) 

A gregarious locust placed to the left of the swarm at a distance equal to the attraction length scale 
a g = 0.14 m should also move to the right, but with a slower velocity which we take to be the minimal 
velocity in a crowd, v™ 111 = 144 m/hr. Thus 

^(-0.14,0) = {-VQ 3 * Pgroup U(x)} | ( _ Q 14 Q) = vf n . (8) 

These two conditions determine R g = 5.13 m 3 / (hr ■ locust) and A g = 13.33 m 3 / (hr ■ locust) In the main 
text, we present numerical simulations of Eqs. ([l} in one spatial dimension. For these simulations, we take 
<5i.2 , r S) Tg, and a g as above, since these parameters do not depend on spatial dimension. For the remaining 
parameters, we follow a process similar to that described above, and choose k\, 2 — k — 8 locusts /m, 
R s = 6.83 m 2 /(hr ■ locust), R g = 6.04 m 2 /(hr ■ locust), and A g = 12.9 m 2 /(hr ■ locust). 



Homogeneous steady states 

For any set of initial conditions, the mean locust density p is known, and corresponds to the total 
density at the homogeneous steady state (HSS). Accordingly, there is a family of homogeneous steady 
states parameterized by pq. The corresponding solitary and gregarious HSS components, obtained by 
setting time and space derivatives to zero in Eqs. ([1]) are 

py5ikl(kl+ pp 
5ik 2 k 2 + 8 Y k\p 2 + 5 2 k 2 p 2 + 5 2 pt 



~ X '-2 ;„2 i X ;„2 „2 i X. U2 n 2 i x „4 ' 



52 P l{k 2 + p 2 ) 

50 5 1 k 2 k 2 + 5 1 k 2 p 2 + S 2 k 2 p 2 + 6 2 pf { ' 

When we later consider stability of homogeneous steady states, it will be convenient to discuss the 
fractions (f> s _ g of solitarious and gregarious locusts, where (j> s + <j) g = 1. Using Eqn. @, we know that for 
homogeneous steady states, 

/ 3o 



■so + .9o 
1 

so/ffo + 1' 



(10a) 
(10b) 



1 + lK 2 ^ \ ■ (10c) 



3 



Here, 7 = 61/82 is the ratio of maximal solitarization rate to maximal gregarization rate, K = fci/fe is 
the ratio of the characteristic solitarization and gregarization densities for individuals, and ip = po/^2 is 
a rescaled density. Note that (j> g is monotonically increasing in ip, and hence in po; that is to say, as total 
density increases, the gregarious fraction increases. 



Linear stability analysis 

To study the stability of the HSS in Eqs. we consider small perturbations S\,g\ about so?5o 

s(x,i) = s + si(x,i), ff(x,i) = g +g 1 (x,t), (11) 

so that p(x, t) = so + 9o + si(x, <) + g\ (x, i). Substituting Eqn. ([TTjl into Eqn. (JTJ and expanding to first 
order in the perturbations, we find the linearized equations 

si = s Q s * V 2 (si +.91) -A Sl +Bg 1: (12a) 

9i = 9oQ g *^ 2 {s 1 +g 1 )+As 1 -Bg 1 , (12b) 

where 

A = / a (po) + /^(po)*o - /i(Po)<?o, (13a) 

S = fM + f'Mto-&{pQW (13b) 

Here, A, B > for all po > since /1 is a monotonically increasing function of po and /2 is a monotonically 
decreasing one. To further analyze the linearized equations, we Fourier expand the perturbations as 

Sl (x,t) = £s q (i)e^ x , fi 2(x,i) = J> q (t)e*»~ (14) 

q q 

We allow for an infinitely large domain so that there are no restrictions on q; in other situations, q must 
be suitably restricted in order to satisfy boundary conditions. Substituting Eqn. (|14jl into Eqn. ([12} yields 
ordinary differential equations for each Fourier mode amplitude. We write these in matrix form, 

^ - L(g) ( S A , (15a) 



dt \GqJ \Gq 

t(„\ = (- s oq 2 Qs(q)-A s q 2 Q s (q) + b\ 
W - \-9oq 2 QM + A -909 2 Q a (9) -Bj- [ ™> 

Here, q = |q| is the perturbation wavenumber, and Q sg (q) are the Fourier transforms of the two dimen- 
sional social interaction potentials, 

«.W = (TT^- < 16 > 

2nR g r 2 g 2irA g a 2 g 

QM = (1+r 2 g 2)3/2 - (l +a 2g2)3/2- ^ 

The eigenvalues Xi_ 2 {q) of L(g) are 

= -9 2 [«oQ.(«)+SbQ«(«)l , A a = -(i4 + B). (18) 
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Since A2 < 0, instability occurs only when Ai > 0. For convenience, we rewrite Ai in terms of the 
gregarious mass fraction <j) g , 



Ai(g) = -po<? 2 (^-<t>g)Qs{q) + 4>gQ g {q) 



(19) 



Now we factor out the attractive part of the gregarious term, namely 



2w A g a* 



(1 



2)3/2 



(20) 



This yields 



M?) 



2 , ^Agd 2 

-poq <p g — 



(l + a 2 q 2 ) 3 / 2 



l-4> g R s rl{l + a 2 g q 2 f/ 2 
^ g A g a 2 g {l + rW)m 



R g r 2 {\ + a 2 g q 2 fl 2 
A g a 2 g (1 + 7-^)3/2 



(21) 



Since the prefactor is negative, and we seek conditions for a positive eigenvalue (signifying growth of 
perturbations, and hence instability), we focus on when the term in square brackets becomes negative. 
The dependence on <fi g occurs via the prefactor (1 — (f> g )/(f> g in front of a positive term. For possible 
instability, this term should be small, meaning that <f> g should be sufficiently large (since this prefactor 
is monotonically decreasing with 4>g)- Since <pg increases monotonically with po (as discussed above), 
instability may occur as po is increased. 

We now show that instability first occurs at the wavenumbcr q = (meaning that perturbations that 
first lead to instability are long wavelength). We again focus on the bracketed quantity in Eq. (12T1) . If 
this term becomes negative, it must do so for the value of q at which the first two terms are (together) 
minimized, since these are positive terms and the negative term, —1, is a constant. It is biologically 
reasonable to assume that a g > r s (with equality achieved for our chosen social interaction parameters). 
Therefore, the first term is either constant or monotonically increasing in q. It is also biologically reason- 
able to assume that a g > r g , in which case the second term is monotonically increasing in q. Thus, the 
first two terms together are monotonically increasing in <j, so their minimum occurs at q = 0, and this 
will be the first wavenumber to trigger instability. Thus, if we are looking for the instability that occurs 
as <p g increases, it is sufficient to consider what happens at q = 0. 



We substitute q = into the bracketed term in Eqn. (f2Tj) and ask for what value of 4> g the resultant 
expression changes sign (to find the threshold level of gregarious locust fraction needed for instability). 
Setting that bracketed term to zero we obtain 



R s ri 



Rgr 2 



A g a 2 g 



(22) 



Instability is achieved for values of <p g greater than this threshold value. 

To obtain a more explicit condition for instability in terms of the density po, we substitute <p* into 
Eq. (|10[) , which relates gregarious fraction to total (scaled) density. Rearranging, we obtain the bi- 
quadratic equation 



where 



A^ + + C = 0, 



K- 1 
— 7if : 



B = 
C = 



-1-7 



(23) 

(24a) 

(24b) 
(24c) 
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For any biologically meaningful solutions, the solution for ip 2 must be positive. From the quadratic 
formula, we have 

Since A > and C < 0, the discriminant is positive. Hence, for the plus sign choice, ^j 2 > 0. For the 
minus sign choice, ip 2 < and hence we eliminate this possibility. The final result for the critical scaled 
density is 



This is the result that we use to produce instability contours in the if -7 plane (Fig. 2 in the main paper). 



Numerical simulation method 

We simulate Eqs. CO)-© in one spatial dimension. We use periodic boundary conditions on a domain of 
length L with a fine grid consisting of N = 1024 points (necessary to resolve the steep edges of clusters 
that form). To approximate an unbounded domain, one may take the limit of large L. The social 
interactions Q s g in ((4]) must be adapted to be commensurate with a periodic domain. We begin with the 
function Q(x) = c^ x ^ r , which is the building block of Q s ,g- We calculate the discrete Fourier transform 
T of —d x Q on our domain as 

T{-d x Q(x)\ = (27 

r cosh(A/r) — cos(Aq) 

where r is the decay length scale in Q and A = L/N is the grid spacing. From Eqn. (|27p it is straight- 
forward to compute the Fourier transforms of Q s . g . Convolutions are equivalent to products in Fourier 
space, providing excellent computational savings (and thus justifying the choice of a periodic domain). 
We compute velocities by convoluting the density with —d x Q s g pseudospectrally. The flux term in 
Eqs. (fT]) is instead evaluated via a fourth-order accurate central finite difference. 

The emergence of discontinuities in s and g causes ringing in the pseudospectral evaluation of the 
velocity term. In order to smooth this effect, we incorporate small amounts of numerical diffusion. 
Another standard approach would be to incorporate high wave number filtering in the simulation. We 
choose numerical diffusion because it also serves as the macroscopic description of random motion, which 
locusts certainly display. We implement diffusion in a split-step manner, alternating with the dynamics of 
Eqs. (HJ-©. Time-stepping is performed with the fourth-order Rungc-Kutta method. We also threshold 
our velocity field at every time step so that it docs not exceed u™ ax . Without this thresholding, individual 
locusts achieve velocities of up to approximately 1.5 times i;™ ax at an intermediate stage of our simulation. 
It is crucial to point out that this thresholding only affects the speed of the transient clumps; it does 
not affect the initial instability (which is small amplitude, and thus has a small velocity) and similarly, 
it does not affect the late-stage bulk dynamics (which are nearly spatially stationary). 
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Abstract 

Locusts exhibit two interconvertible behavioral phases, solitarious and gregarious. While solitarious indi- 
viduals are repelled from other locusts, gregarious insects are attracted to conspccifics and can form large 
aggregations such as marching hopper bands. Numerous biological experiments at the individual level 
have shown how crowding biases conversion towards the gregarious form. To understand the formation of 
marching locust hopper bands, we study phase change at the collective level, and in a quantitative frame- 
work. Specifically, we construct a partial integrodiffcrcntial equation model incorporating the interplay 
between phase change and spatial movement at the individual level in order to predict the dynamics of 
hopper band formation at the population level. Stability analysis of our model reveals conditions for an 
outbreak, characterized by a large scale transition to the gregarious phase. A model reduction enables 
quantification of the temporal dynamics of each phase, of the proportion of the population that will even- 
tually gregarize, and of the time scale for this to occur. Numerical simulations provide descriptions of 
the aggregation's structure and reveal transiently traveling clumps of gregarious insects. Our predictions 
of aggregation and mass gregarization suggest several possible future biological experiments. 

Author Summary 

Locusts such as Schistocerca gregaria, Locusta migratoria, and Chortoceites terminifera periodically form 
highly destructive plagues responsible for billions of dollars in crop losses in Africa, the Middle East, 
Asia, and Australia. These locusts usually exist in the so-called solitarious behavioral phase and seek 
isolation; gregarious individuals, however, are attracted to conspecifics. Previous experimental work has 
uncovered the causes of phase change in individual insects: principally, sustained exposure to sparse or 
crowded conditions. An open problem is to understand the intrinsic roles that phase change and social 
interaction play in the transition from an initially disperse, solitarious population to an aggregated, 
destructive marching hopper band of gregarious individuals. To this end, we construct a mathematical 
model that describes the interplay of phase change and spatial dynamics. Through analysis and numerical 
simulations, we determine a critical density threshold for gregarious band formation and quantify the 
collective phase change over time. We also discuss implications of our work for preventative management 
strategies and for possible future biological experiments. 

Introduction 

Outbreaks of locusts such as Schistocerca gregaria, Locusta migratoria, and Chortoceites terminifera 
regularly afflict vast areas of Northern Africa, the Middle East, Asia, and Australia. Depending on 
climate and vegetation conditions, billions of voracious locusts aggregate into destructive swarms that 



2 



span areas up to a thousand square kilometers. A flying locust swarm can travel a few hundred kilometers 
per day, stripping most of the vegetation in its path [H3]- A recent locust plague in West Africa (2003- 
2005) severely disrupted agriculture, destroying $2.5 billion in crops destined for both subsistence and 
export. Despite control efforts totalling $400 million, loss rates exceeded 50% in certain regions [SJ|S]. 
These numbers alone attest to the urgency of finding better ways to predict, manage, and control locust 
outbreaks. 

Between outbreaks, locusts are mainly antisocial creatures who live in arid regions, laying eggs in 
breeding grounds lush with vegetation. Resource abundance may, on occasion, support numerous hatch- 
ings, leading to a high population density. Overcrowding at resource sites promotes transition to a social 
state in a self-reinforcing process. The social locust nymphs may display mass migration behavior. Within 
the newly formed group, individuals cohere via sensory communication, whether visual, chemical, and/or 
mechanical [3], Outbreaks may be exacerbated in periods of drought, when large numbers of locusts 
congregate on the same breeding or feeding grounds [7HS]. 

Locusts are phase polyphenic: while sharing the same genotype, individuals may display different phe- 
notypes plIIfTT] that incorporate variations in morphology [12] , coloration [13] , reproductive features [T4] 
and, significantly, behavior [151116] . An individual can change from a solitarious state (preferring isolation) 
to a gregarious one (seeking conspecifics). Behavioral state is plastic [3]lll [ [l"5" ] and strongly dependent 
on local population density: in sparse surroundings, a gregarious locust transitions to the solitarious 
state [15] and vice versa in crowded environments. These phase transitions are called solitarization and 
gregarization. Gregarization dominates when large numbers of locusts gather at the same site, potentially 
leading to a destructive outbreak [8j[9]. 

Locust gregarization may be induced by visual, olfactory, or tactile cues. For the desert locust 
Schistocerca gregaria, the most potent stimulus is tactile: repetitive stroking of the femora of hind 
legs [15Tri7] functions as a crowding indicator. Mechanosensory stimulation of leg nerves leads to serotonin 
cascades in the metathoracic ganglion, and initiates gregarious behavior [TBTfTS] . Gregarization can be 
induced by rubbing a locust's hind leg for 5 s per minute during a period of 4 hr [17] . Cessation of 
physical contact leads to solitarization after 4 hr, though the degree of solitarization achieved during 
that time depends on the individual's ancestry. 

Experiments and models have shed much light on how group alignment [19H22] and group motion 
[231124] depend on group size or density and treatments such as diet and denervation. For instance, a low- 
protein diet (which motivates cannibalism in locusts) leads to stronger interactions between individuals 
and lowers the threshold density beyond which mean speed and group coherence increase [24] . Other 
data-driven studies include models based on a well-known physics paradigm for self-propelled particles [25] 
and explore the transition between a disordered and a coherent marching group. Both [26] and [27] study 
the dynamics of rolling patterns formed by flying, gregarious swarms. A logistic map was introduced 
in |28] to describe phase change via a birth rate and a carrying capacity dependent on population density 
modulated by stochastic effects. 

Our current work complements previous locust modeling studies in several ways. First, many of the 
previous models are individual-based (Lagrangian) simulations, where the position, velocity, and interac- 
tions of individual locusts are tracked [191 20 ,22 24] . Ours is density-based (Eulerian), allowing techniques 
of partial differential equations (PDEs) and their extensions (integro-PDEs) to be utilized. Second, we 
concentrate on gregarious-solitarious transitions not yet explicitly considered in j2Qi[24] . We address in- 
trinsic attractive-repulsive social interactions, whereas many current models consider interactions with 
clumped resources and environmental heterogeneity as their focal points [HE] . Finally, some models [24] 
include anisotropic interactions such as different responses to anterior and posterior neighbors, or consider 
Newtonian dynamics. To explore minimal mechanisms sufficient for band formation, our work instead 
uses isotropic interactions and a kinematic approach. The open problem we address via mathematical 
modeling is to quantify and describe collective gregarization, a key, early process that necessarily occurs 
before the emergence of a destructive locust outbreak. We do this by linking the physiology of individual- 
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level phase change and interphase interactions to predictions at the level of the gregarious hopper band 
as a whole. 

We investigate the onset of an outbreak by constructing a continuum mathematical model of behav- 
ioral phase for interacting gregarious and solitarious locusts. We classify and quantify group dynamics 
in wide swaths of parameter space, a task which is challenging by numerical techniques alone. We find 
that in the limit of low densities, both phases arc uniformly spread and the solitarious phase dominates. 
For sufficiently large populations, a dense, traveling patch of gregarious locusts suddenly emerges, while 
solitarious locusts become more and more scarce. We identify locust clustering at high densities with the 
onset of a hopper band. Through analysis of our model, we calculate the critical density beyond which 
the gregarious group forms, and for the final ratio of gregarious to solitarious locusts. We determine these 
quantities in terms of behavioral parameters at the level of individual locusts, hence connecting individ- 
ual and group properties. Our model also displays population-level hysteresis, which has implications for 
locust management. 



Model construction 

Locusts in a group are subject to attractive and/or repulsive forces based on combined sensory, chemical, 
and mechanical cues that affect their motion. We assume that sensing is directionally isotropic, a reason- 
able approximation ;29i for organisms receiving sensory inputs of a variety of types, although directional 
models are possible as well [3D]. Rather than tracking individual locusts, we consider a population den- 
sity field p(x, t) moving at velocity v(x, t). Continuum population modeling |31H32| allows us to apply 
analytical tools in order to characterize swarm formation and structure. Our work draws from classic 
swarm modeling in which a conserved population density field p moves at a velocity v that arises from 
social interactions: 



This is the well-known mass balance equation that tracks individuals moving collectively at velocity v. It 
is typically assumed that individuals can sense the population density nearby, and that this sensing gives 
rise to attractive-repulsive social forces F, or alternatively, social potentials Q (the negative gradients of 
which are forces). Within this context, the contribution p(x , t) of a small clump of individuals at location 
x' to the force on the individual at position x is given by F(x — x')p(x/ t) = — VQ(x — x')p(x', t). The 
corresponding velocity is proportional to the forces exerted by neighbors at all spatial locations, so that 
v(x, t) is given by integration over all x' as 



The expression for the velocity v(x, t) in Eqn. ^ is a convolution of the density p(x, t) and the social 
interaction force — VQ(x — x'), which describes the influence of the locust population at location x' 
on that at location x. This is a common formulation of so-called nonlocal interaction models |33H36) . 
which capture interactions that are spatially distributed, in contrast to pure partial differential equations, 
which include only local terms such as derivatives and gradients, and which describe interactions only 
over infinitesimal ranges. Nonlocal aggregation models have been studied for various social interactions 
Q; known solutions include steady swarms, spreading populations, and contracting groups. We use the 
notation v = — VQ * p to denote the convolution in Eqn. ©. We assume that Q(x — x') is radially 
symmetric and depends only on the distance between x and x'. The detailed forms of Q in the case of 
solitarious and gregarious locusts will be described later. 

To adapt Eqs. ([1} and (0) to biphasic insects, we introduce separate density fields for solitarious and 
gregarious locusts, s(x, t) and <?(x, i), respectively, and the total local density p = s + g. With marching 
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p t + V • (pv) = 0. 



(1) 




(2) 
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locusts in mind, we consider a two-dimensional geometry, with 51 representing the spatial domain and 
x = (x, y) as spatial coordinates. We now include the phase transitions between solitarious and gregarious 
locusts. To do so, we define two density-dependent functions, fx(p) for the the rate of gregarious-to- 
solitarious transition, and f% (p) for the rate of solitarious-to-gregarious transition. Our model thus reads 



s + V • (v.«)= -h{p)s + h(p)g, (3a) 

.9 + V • (v 3 . 9 )= h{p)s- h{ P )9, (3b) 

where the velocities are given by 

v s = -W(Q s *p), w g = -V{Q g *p). (4) 

These equations are complete once we specify the solitarious and gregarious social interactions Q s ,g and 
the density-dependent conversion rates /i,2- Since solitarious locusts are crowd-avoiding, we take Q s to be 
purely repulsive. Gregarious locusts, on the other hand, are attracted to others, except for short-distance 
repulsion due to excluded volume effects. Hence, we model Q s and Q g as 

Q s (x-x') = i? s e-l x - x 'l/^, Q 9 (x-x') =i? 9 e-l x - x 'l^ -A 9 c-' x - x 'l/ Q 9 , (5) 

where R s , R g , A g are interaction amplitudes that determine the strengths of attraction and repulsion, 
and r s , r g and a g are interaction length scales that represent typical distances over which one locust can 
sense and respond to another. 

The above forms of Q s ,g describe social interactions that decay exponentially away with distance 
from the sensing individual and are chosen to be isotropic for simplicity. As evident from Eqn. ([5]), Q s is 
purely repulsive for all choices of R s and r s . On the other hand, Q g is the difference of two exponentials, 
implying that there may be a distance at which repulsion and attraction balance, resulting in no net 
contribution to the velocity. The location of this balance point can be obtained by imposing —\7Q g = 
to obtain the critical distance 

d= _Vi_] n f5!5ay (6 ) 

a g - r 9 V A a r 9 J 

Depending on the choice of social interaction parameters, the expression for d may yield unphysical 
results such as negative distances. The distance d also pertains only to two isolated locations x and x' 
and does not capture population-level features. Even for meaningful values of d, a collection of individuals 
interacting under Q g may disperse, aggregate, or clump. It is thus important to choose the appropriate 
parameter ranges for a g , r g , A g and R g so that the tendency of gregarious locusts to aggregate is modeled 
properly. Mathematical studies have shown that in order for cohesiveness to occur, the parameters in 
Q g must lie in a particular regime that leads to clumping J37]. Thus, we require R g a g — A g r g > so 
that repulsion dominates at short length scales, and A g a 2 g — R g r g > so that attraction dominates at 
longer ones. Taken together, these conditions guarantee a meaningful critical distance d and macroscopic 
clumping behavior. We assume these conditions to hold for the remainder of this paper. 

It remains to specify how density affects transitions from one phase to another. We call upon the 
biological observation that at higher densities, gregarization proceeds more quickly and solitarization 
more slowly. We model the phase conversion rates with the rational functions 

/l(p) = m= W*») a (7) 

l + (p/fcl) 2 l + (p/fc 2 ) 2 



The parameters 5\ t 2 are maximal phase transition rates and are characteristic locust densities at which 
fi^2 take on half of their maximal values. Note that fi decreases with p, capturing the inverse relationship 
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between solitarization rate and density, while /2 increases with p and saturates at 82 , describing speedier 
gregarization at higher densities. 

Our complete model consists of Eqs. ©-(CO) together with initial conditions specifying s(x, 0) and 
<?(x, 0). We consider a spatially periodic domain, which simplifies both numerical simulation and math- 
ematical analysis. In certain laboratory studies using ring-shaped arenas, such boundaries are natural 
(while being less ideal for comparison with field studies) [20j . We do not include locust reproduction or 
death as these occur on much longer time scales than phase change. 

The model presented here is a general one containing some fundamental elements of locust dynam- 
ics. This work can be readily modified and extended to include details pertaining to different locust 
species, interactions with the surrounding environment, locust reproduction, and more. For instance, 
in our model, we have not explicitly accounted for the differing activity levels of solitarious and gre- 
garious individuals [11) . Additionally, while gregarization is relatively fast for Schistocerca gregaria, full 
solitarization may occur only after several generations of locusts. The phase conversions of Chortoicetes 
terminifera, on the other hand, are characterized by similar timcscales for the two phase conversions, 
so that both gregarization and solitarization occur rapidly within the lifetime of a single locust individ- 
ual [38]. On another note, vegetation or waterway patterns may impose spatial inhomogeneities such as 
non-uniform initial distributions of solitarious locusts, or attraction to preferred sites. Preexisting models 
in the literature have pointed out the important link between the spatial distribution of vegetation, as 
well as nutritional quality, on locust clustering, gregarization, and swarming [SHU [Ml SO]- All of these 
elements could be used to refine our model for predictive purposes. However, as the first work in the con- 
tinuum modeling of locust population phase change, ours begins with the fundamental model contained 
in Eqs. ([3])- (J7J). Our model is complementary to the preexisting ones in that we focus on how inherent 
inter-individual interactions can lead to gregarization and swarming, even in a spatially homogeneous 
environment. Multi-generational dynamics, differential activity levels, resource distribution, and related 
factors could be considered as possible extensions of our work. 

Parameter selection 

Some of our results are analytical formulas, which may be evaluated for any desired parameter values. 
Other results depend on numerical computations, and these require specific choices of parameters. For 
these results, we consider two different sets of phase transition parameters. (1) Most of our numerical 
results have been obtained using our default set of parameters, based on estimates from the biological 
literature. Specifically we take 81,2 = S = 0.25 hr^ 1 , corresponding to a gregarization time scale of 
approximately 1/8 = 4 /i?' for desert locusts (for whom some - but not total - solitarization occurs on 
the same time scale) [TTlll7) . We also take k\^ = k = 65 locusts / m 2 , since for desert locusts, the critical 
density for the onset of collective motion is 50 - 80 locusts /m 2 [24]. We will allow for some deviation 
from Si ~ 82 and k\ = ^2 via a parameter sensitivity analysis. (2) To examine situations with large 
differences in the rates of gregarization and solitarization, we consider an alternative set of parameters 
with Si = 0.025 hr~ l and 82 = 0.25 for -1 , so that gregarization is an order of magnitude faster that 
solitarization. We take k\ = 20 locusts / m 2 and &2 = 65 locusts / m 2 to model a gregarious-to-solitarious 
transition that occurs at a higher density threshold than the solitarious-to-gregarious transition. 

We use the same social interaction parameters for all results (variations from this set are accounted for 
by a sensitivity analysis). To estimate the social interaction length scale parameters in Eqs. ([5]), we apply 
the results of [3DK21], which identify the "sensing range" of a desert locust as 0.14 m, and the "repulsion 
range" as 0.04 m, close to the approximately 0.05 m body length of a desert locust at the fifth instar 
of its development. For the gregarious phase we thus set the repulsion length scale at r g = 0.04 m and 
the attractive one at a g = 0.14 m, corresponding to the experimental sensing range. These choices agree 
with theoretical studies showing that for cohesive swarms, attraction occurs over longer length scales 
than repulsion [33][3T|. We also assume that solitarious locusts are repelled from others at their sensing 
range, so that r s = 0.14 m. These choices satisfy r g < a g = r s which is assumed for the remainder of 
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this paper. 

Finally, we estimate R s , R g , and A g via explicit velocity computations. The speed of a locust when 
it is alone varies between 72 - 216 m/hr, depending on diet |24| . At the upper end, this is roughly one 
body length per second. When it is moving in a group, the individual's speed varies in a tighter range of 
144 - 216 m/hr [24]. In making our phase-dependent velocity estimates, we interpreted the "moving alone" 
and "moving in a group" data as typical to solitarious and gregarious locusts, respectively. Using these 
biological measurements and Eqn. (j4j, we find R s = 11.87 m 3 /(hr ■ locust), R g = 5.13 m 3 /(hr ■ locust), 
and A g = 13.33 m 3 /(hr- locust). Details are given in Text SI. Our choices of social interaction parameters 
satisfy conditions mentioned in the previous section, namely R g a g — A g r g > 0, and A g a 2 g — R g r 2 > so 
that gregarious insects will clump. 

Most of our parameter choices have been inferred or estimated from published laboratory experiments. 
It is possible however, that in the field, some parameter values may be quite different from the ones we 
have used. For instance, locusts in the field may pause while marching to perch on the vegetation, giving 
rise to an effective speed that is lower than what measured in lab experiments, where perching docs not 
occur. It is also noteworthy that gregarious locusts are more active than solitarious locusts, a fact that 
is reflected by our method of choosing R s ,R g ,A g from estimates of the velocities of individuals when 
moving alone and in a group. As we describe below, we analyze our model varying all parameters within 
reasonable bounds: our results are qualitatively the same. 

Results 

We first determine the simplest solutions to the model, namely those for which the densities of gregarious 
and solitarious locusts are in a spatially uniform steady state. We probe the stability of that uniform state 
using linear stability analysis (LSA), a calculation that addresses whether small, spatially nonuniform 
perturbations grow or decay. This is equivalent to determining the signs of eigenvalues of the linearized 
system, where positive (negative) eigenvalues imply growing (decaying) perturbations. The rate of ini- 
tial growth/decay depends on the wavenumber of the perturbation. The growing perturbations can be 
interpreted in terms of nascent aggregates of locusts, and the wave numbers as the number of aggregates 
per unit area. The analysis provides a condition for the onset of aggregation, namely the emergence 
of positive eigenvalues of the linearized model. In our case, this aggregation condition is shown below 
in Eqn. (|14[) . LSA cannot, in general, predict the ensuing dynamics once perturbations have grown to 
a large size. Further analysis uses an approximation to eliminate the spatial dependence of the model, 
which enables an analytical prediction of the proportion of solitarious and gregarious locusts on a longer 
time scale. To visualize the dynamics of aggregation, we perform numerical simulations in one spatial di- 
mension using the linear stability analysis to identify regimes of interesting behavior. The model displays 
population-level hysteresis. 

Homogeneous steady states 

The solitarious sq and gregarious go homogeneous steady-state (HSS) solutions of Eqn. ([3]) can be written 
in terms of the total uniform density po, which is simply the mean value of p for a specified initial 
condition. The full expressions for so and go in terms of po appear in Text SI; in the small po limit these 
are approximately 

^23 ^23 

So * Po - ^f ' 9o ~sD4 Po > (8) 

while in the limit of large po we find 
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The low density HSS is thus composed mostly of solitarious locusts and vice versa for the high-density 
case, showing the non-monotonicity of sq with respect to total density po- In Fig.QTA) we plot the HSS 
so (middle solid blue curve) and go (middle broken green curve) for our default set of phase change 
parameters, k = 65 locusts / m 2 and 8 = 0.25 hr" 1 . 

As shown, so initially increases with po- At a critical density p*, so reaches a maximum, whereas 
go keeps increasing monotonically. Fig.[TjB) shows a blow-up of the region near p*. For our default 
parameters, the maximum value s™ 3 '* is attained at p* = k, the same density value for which solitarious 
and gregarious densities coincide so that s nax = so(p*) = go(p*) = k/2. However, this feature is a result 
of our choice k\ = ki and 6\ = 62- In general, the point of maximum solitarious density and the point of 
equal solitarious and gregarious density do not coincide, as is directly deducible from the full expressions 
for so and go in Text SI. To give a sense of detuning from our parameter estimates, we also calculate and 
plot so and go for parameter sets chosen randomly from uniform distributions centered at our estimated 
default set of values for {81, 62, ki, k%}. The bottom and top curve in each set show the 25th and 75th 
percentile values. 

We also study a much more general case where <5i ^ ^2,fci ^ fe, in keeping with the distinct rates 
of transition and critical transition densities seen biologically. As an alternative way to understand the 
HSS solutions, we consider the fractions (p s>g of solitarious and gregarious locusts, where <p s + 4> g = 1. As 
shown in Text SI, for the HSS, 



Here, 7 = 81/62 is the ratio of maximal solitarization rate to maximal gregarization rate, K = fci/fe is 
the ratio of the characteristic solitarization and gregarization densities for individuals, and ip = po/^2 is 
a rescaled spatially homogeneous density. The gregarious fraction <f> g is monotonically increasing in if), 
and hence in po; that is to say, as total density increases, the gregarious fraction increases. For small po, 
4> g « 0, but as po increases, there is a crossover between solitarious and gregarious populations. Uniformly 
spread solitarious populations cannot be sustained when the density is too high: the gregarious state will 
necessarily become the dominant one. 

Linear stability analysis 

To determine conditions under which a nearly uniformly spread locust population aggregates or disperses, 
we study the linear stability of the HSS (details appear in Text SI). The calculation is a standard but 
somewhat tedious exercise. In nonlocal systems such as ours, linear stability results depend on the 
Fourier transforms Q s ^ g (q) of the interaction potentials Q s , g . For our locust model, the stability of the 
HSS depends on the eigenvalue 



where q = |q| is the perturbation wave number and the Fourier transforms Q Sig (q) in two dimensions are 




(10) 



= -q 2 soQs{q) + goQg(q) 



(11) 



Qs(q) 



2irR s r 2 s 



(12) 



Q 9 (q) 




(13) 



( 1+r 2 9 2)3/2 (l + a 2 g 2)3/2- 



Observe that the eigenvalue Ai (q) depends on all of the individual-based parameters governing rates of 
phase change (via sq and go) and all of the social interaction amplitudes and length sensing length scales. 
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The HSS derived in the previous section is stable to small perturbations if \i(q) < for all q. If Ai (q) > 
for some q, then the HSS is unstable to perturbations of those wave numbers. 

Our full analysis of this eigenvalue appears in Text SI. We formulate the instability condition in 
terms of <fi g , 

If this condition is satisfied, initially small perturbations from the uniform steady state will grow. This 
inequality is a key result, and implies that if a sufficiently large fraction of the population is gregarious, 
the HSS solution is unstable. To obtain a more explicit condition in terms of the density po, one must 
substitute <f>* into Eqn. (flU)) . which relates gregarious fraction to total (scaled) density. One may then 
calculate the critical density p above which the HSS is unstable. Since 4> g and po are monotonically 
related, we conclude that the HSS solution is unstable for sufficiently dense populations. The algebra is 
tedious, and relegated to Text SI. Instead, we present a contour plot in Fig.[2]which succinctly illustrates 
the stability features of the HSS. The phase change parameter ratios 7 = 81/82 and K = k\/ki vary 
along the horizontal and vertical axes and the contours indicate the critical value of rescaled density 
tp* = For scaled densities greater than ip*, the HSS solution is unstable. The critical scaled 

density is monotonically increasing in both 7 and K. (Note that for an accurate biological interpretation, 
one must multiply ijj* by &2 in order to obtain the unsealed critical density p$.) 

Upon inserting our default parameters in Eqn. (|14[) we find that the homogeneous solution is unstable 
for po > Pq = 62.3 locusts/m 2 . This value corresponds to the left border of the grey region in Fig.[T] For 
Po > Pq, to the right of the border, we expect the onset of a locust hopper band, i.e., formation of patches 
of high locust density that can seed the clustering and gregarization of other locusts. In Fig.[IJ linear 
instability can occur even at densities po for which sq exceeds go for our chosen parameters (represented 
by the center solid blue and center broken green curves). This result implies that the onset of instability 
leading to mass gregarization can take place even if solitarious locusts initially outnumber gregarious 
ones. We will later discuss mass gregarization in more detail. To visualize detuning from this set of 
parameters, we include the 25th and 75th percentile values of Po for onset of instability as vertical purple 
lines; these are again calculated by drawing 10,000 random samples of the parameters k\^, <5i,2, Rs.g, 
r s g , A g , and a g . As seen from Fig.fljb) our conclusions are robust across the randomly chosen parameter 
sets. 

For our default set of biological parameters, <fi* w 0.479 via Eqn. ([M)) and Po turns out to be near 
k = fci.2. We stress that generically, it is not the case that Po needs to be near fci and/or k-i- For our 
default parameter set, K = 7 = 1, in which case ip = 0.959, so that the critical value Po is 95.9% of 
^2, namely 62.3 locusts/m 2 . However, for different choices of K and 7, drastically different outcomes 
are possible. For instance, for our alternative parameter set where K ps 1/3 and 7 = 1/10, the critical 
density is po = 15.9 locusts/m 2 , which is quite disparate from the individual gregarization density 
of 65 locusts/m 2 , and is also less than the solitarization density of 20 locusts/m 2 . Furthermore, for 
different choices of the social interaction parameters entering into Eqn. (|14|). it is possible to obtain a 
critical gregarious fraction </>* that is much less than 1/2, meaning that instability and clumping can 
occur even with just a few gregarious insects. 

For po > Po, we can also find the wave number q max corresponding to the most rapidly growing 
perturbation. Fig. |3] shows g max for our chosen parameters (center curve) as well as the 25th and 75th 
percentile values over the 10,000 random parameter draws. The most unstable wave number g max grows 
rapidly as a function of po and then saturates at q max ~ 8.89 m _1 , corresponding to a length scale 
27r/<2 , max ~ 0.71 m and indicating that the most quickly growing perturbations occur on the length scale 
of a few locust bodies. 

Our linear stability analysis describes the behavior of small perturbations of uniform steady states, 
and is not expected to predict long-term or large-amplitude dynamics. For large perturbations, linear 
analysis is void. Additionally, even to analyze small perturbations of states other than uniform steady 
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states, a different analysis would be needed. 
Numerical simulation 

To illustrate the swarm dynamics described by Eqn. (|3|), we simulate the model on a one-dimensional 
periodic domain of length L = 3 m for a total population of M = 50 locusts. Periodicity of the domain 
is an important aspect of a robust numerical platform devised for these simulations: we exploit the 
fact that convolutions Q * p are easy to compute in Fourier space (where they are simply products, i.e., 
Q-p), which significantly reduces the computational overhead. Computational issues associated with such 
convolutions also restrict us to one-dimensional simulations at present. At t = all locusts are solitarious 
and are randomly perturbed from the uniform density s = M/L, where M is the total population mass 

M= I pdx. (15) 
Jn 

We adjust some parameters so as to adapt our model to the one-dimensional case. Specifically, one must 
take square roots of k\^ in order to collapse densities in a square to densities along a line segment. 
Consequently, for our default parameter set we choose fci j = k = 8 locusts /m and 6i t 2 = S = 0.25 hr~ x , 
whereas for the alternative set we use k\ = 4.5 locusts/m, &2 = 8 locusts/m, Si = 0.025 hr^ 1 and 
62 = 0.25 hr -1 . In both cases we take the interaction amplitudes R s = 6.83 m 2 /(hr ■ locust), R g = 
6.04 m 2 /{hr ■ locust), and A g — 12.9 m 2 /(hr ■ locust), which have also been adapted from their original 
values to the one-dimensional case. The interaction length scales r s , r g , and a g are the same as for the 
two-dimensional case. Details of the numerical method and the parameter choices appear in Text SI. 

Results are shown in Fig. |4] for the default parameter set and in Fig. [5] for the alternative set. In each 
case, the snapshots show s(x,t) (dashed blue curve) and g{x,t) (solid green curve) at selected times. 
Starting from the randomized solitarious state at t = hr, locusts rapidly redistribute to a roughly 
spatially uniform density until t w 3 hr. Tiny variations are present but not visible on the scales of these 
figures. Grcgarization and subsequent rapid spatial segregation follow. In Fig. [4) between t « 3.42 hr and 
t w 3.47 hr, two compactly supported clumps of gregarious locusts emerge, superposed on a background 
of sparse, solitarious individuals. A similar transition occurs between t m 3.15 hr and t w 3.17 hr 
in Fig.[5j but for these parameter values, we find initial clustering with three, rather than two density 
peaks. The number (or alternatively, length scale) of transient clumps that form appears to be selected 
dynamically. This intermediate dynamical selection process and the coarsening that ensues are avenues 
for future numerical and analytical investigation. In each example, the disjoint clusters quickly merge due 
to the long-range attraction of gregarious individuals. A single remaining pulse is formed by t « 3.49 hr 
in both cases and travels until t « 6.5 hr, at which time the majority, but not all, of the solitarious locusts 
have transitioned to the gregarious form. Grcgarization continues during the subsequent hours, albeit 
at a slower rate. For both figures, the grcgarization of the final clump continues slowly, approaching an 
equilibrium at exponentially long times. 

To study the locust grcgarization process further, we define the total mass of solitarious and gregarious 
locusts, S and G, as 

S= I sdx, G = I gdx, (16) 
Jn Jn 

so that the total population mass is M — S + G. Wc also define the mass fractions 

cb s = S/M, cb g = G/M, cp s +4> g = l, (17) 

which wc before calculated for HSS solutions, but we now generalize for spatially varying states. These 
quantities will be useful to further our mathematical analysis. Fig. [6] shows 4> s (t) (blue curve) and <f> g (t) 
(green curve) as arising from the numerical simulations depicted in Fig. [4] and Fig.[5l Several distinct 
regimes are visible, and we discuss these below. 
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Spatially-homogeneous and spatially-segregated bulk theories 

As visible in the second and third panels of Fig. 2] and Fig.0 the early-time dynamics of Eqs. are 
approximately spatially homogeneous. As a result, spatially-dependent terms in Eqs. ([3]) are negligible, 
p is approximately constant, and hence the governing equations are linear ordinary differential equations 
(ODEs) that are easily solved. We write the solution of these ODEs in terms of the mass fractions <j) s g , 

= f ( TX7 ( \ I 1 " e- [Mp ° )+Mpo)]t \ , Mt) = 1 - Mt), (18) 
h\Po) + h\Po) 1 } 

where we have used the initial condition <f> s {t = 0) = 1. This analytical solution is plotted in Fig. [5] as a 
dotted line, and agrees closely with the numerical results for the first few hours. 

In the later panels of Fig. [4] and Fig.[5l gregarious and solitarious locusts spatially segregate into areas 
with disjoint support. This means that in each distinct region, p(x,t) ~ s(x,t) or p(x,t) rj g(x,i). We 
thus consider a bulk model reduction to study the dynamics of the two non-overlapping solitarious and 
gregarious populations. In particular, we assume that solitarious locusts are spread throughout most of 
the domain O, covering an area denoted a Sl whereas gregarious locusts are confined to a region with area 
a g . Within these areas, local densities are approximately S/a s and g = G/a g . By integrating Eqs. Q 
over the domain and assuming that s and g are approximately constant in their support, we obtain 



dt 1 + c 2 (j) 2 s 1 + ca4> 2 dt 



(19) 



where 



S 2 M 2 M 2 M 2 

Cl = ^kj> C2 = ^r C3 = <5l; C4 = ^Fr (20) 

The numerical solution of these ODEs (dashed lines in Fig.|6]) agrees closely with the late time full- 
scale numerical simulation results, where we use values of a St9 measured empirically from the terminal 
equilibrium. One can reduce Eqn. (|19[) to a single nonlinear ODE using <f> s = l — <f> g , though this equation 
is not amenable to analytical solution. Since we are interested in the large population limit for which we 
expect potential large scale gregarization, we instead study Eqs. (|19|) for large M. In this case, to leading 
order in M, the bulk model reduces to 

<t> s = S 2 4>s H — = -0 S - (21) 

C4<Pg 

Given the expressions for 03,4 and the fact that ^ S) 0g < 1, the first term is 0(1) whereas the second 
one is much smaller, 0(1/ M 2 ). For large M then, and to leading order, 4> s decays exponentially in time 
with rate 5 2 . This result is based on the assumption of a segregated state, and thus would be expected 
to occur only once segregation is nearly complete. 

Since for large M (nearly) the entire population will eventually become gregarious, the critical density 
Pq is a crucial result. If the population is in the stable regime (where po < Pq) then mass gregarization can 
be avoided and solitarious and gregarious locusts can coexist as uniformly spread populations. However, 
as soon as the population shifts beyond the border of stability (where po > p$) the group grcgarizes and 
the onset of a locust hopper band is inevitable. 



Phase change and hysteresis 

The biological literature discusses the importance of hysteresis in locust phase change, as reviewed, for 
instance, in [11]. It is important to disambiguate the possible meanings and interpretations of phase 
change hysteresis, to place this phenomenon within the context of our model, and most especially, to 
distinguish between hysteretic features at the individual and population levels. 
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One type of hysteresis is simply defined as "rates of gregarization [that] differ from rates of solita- 
rization" [11] . Within our model, this type of hysteresis may be interpreted as cases where Si ^ 82 
or ki 7^ fe- Our results thus far have accounted for this type of hysteresis in three ways. First, for 
our primary parameter set in which 81 = 82 and k\ = k2, we have allowed deviations from equality by 
performing a sensitivity analysis incorporating variations of up to 30% from the base parameter values, 
as represented in the results of Fig.Q] and Fig.0 Second, for our alternative parameter set, we have 
chosen 82 = 10<5i and k2 ~ 3fci. And finally, for analytical results such as the homogeneous steady states 
and their stability, we have obtained analytical formulas into which any values of 81,2 and fci ; 2 can be 
substituted. 

Another interpretation of hysteresis relates to "solitarization [having] two phases: an initial rapid 
phase and a second, slower phase that requires insects to be maintained in isolation across successive 
moults - or generations" |llj . Our model is constructed on the time-scale of a single generation, and 
thus we cannot account for this type of hysteresis, which would require a multi-generational model. 

Finally, we can consider population-level hysteresis. In the context of our model, this type of hysteresis 
refers to macroscopic properties of solutions of Eqn. (J3j> (which are outputs of the model) as opposed to 
differences in individual-level parameters (which are inputs to the model) as in the first type of hysteresis 
described above. Numerical results suggest that our model has population-level hysteresis; see Fig. [7] 
This figure shows the gregarious mass fraction <f> g as the average density po (total mass M divided by 
domain length L) is varied as a control parameter. All phase change, social interaction, and physical 
domain parameters are the same as in Fig. [5] 

The solid (dashed) red curve is an analytical result, representing the stable (unstable) HSS solution, 
as calculated previously via linear stability analysis. For small values of po, the HSS is stable to small 
perturbations. If locusts join the initially stable population the average density p a will increase (assuming 
a fixed spatial domain), shifting the uniform state to the right along the red curve; as yet no clustering 
will be evident. Beyond the point labeled with an asterisk, the uniform HSS loses stability and clustering 
occurs, as previously described. This corresponds to a jump represented by the vertical black arrow. The 
clustered state (green) is now stable. We next ask what happens if locusts are now removed from the 
aggregate, which corresponds to a reduction in po (moving to the left in Fig-tZj - We answer this question 
numerically, by gradually subtracting mass from the population, allowing the system dynamics to evolve, 
and plotting the gregarious fraction as a function of mass. As the mass is slowly removed, the solution 
tracks leftwards along the green curve, indicating the persistence of the gregarious band. In fact, the 
band persists even partway into the regime where the HSS is linearly stable. 

This dynamically observed hysteresis suggest that (for our model) a gregarious aggregation cannot be 
eliminated by reducing overall density to a low enough level where the HSS is linearly stable. This result 
has implications for locust control, as we discuss below. 

Discussion 

In this paper, we derived, analyzed, and simulated a model for the movement, social interactions, and 
density-dependent interconversions of the solitarious and gregarious forms of phase polyphenic locusts. 
The model is based on experimental observations and measurements, parameter values inferred from pre- 
existing work, and basic assumptions about individuals' rules of behavior. We included social exchanges 
via repulsive and/or attractive interactions for gregarious and solitarious individuals, and we accounted 
for phase change with density-dependent transitions, with crowding favoring solitarious-to-gregarious 
conversions. Our model was formulated in terms of continuum equations, allowing us to apply classical 
techniques such as linear stability analysis and bulk approximation. Since these methods were applied in 
two spatial dimensions, our results are relevant to insects aggregating in two dimensional structures such 
as hopper bands. We also provided example simulations in one spatial dimension as proof of principle, 
and as an indication of typical dynamics. 
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Our model explicitly takes into account intrinsic social interactions between individuals, in contrast to 
pre-existing models that focus on how insects respond to quality and spatial heterogeneity of nutrition or 
other environmental factors 8 , 9 , 24 . These approaches are complementary, showing that both intrinsic 
and extrinsic factors that affect local densities also affect the gregarization transition. 

Many of our results are achieved via mathematical analysis. The power of mathematical analysis is 
that it creates an explicit connection between individual-level and group-level quantities, e.g., via the 
inequality Eqn. (|14[) . Once we identify the sensing range and interaction strength parameters in Eqn. (JSJ) 
which govern individual locust attraction to and repulsion from others, we are able to calculate the critical 
density beyond which mass gregarization occurs. 

Briefly, our results and predictions can be summarized as follows: (1) Locusts exist in a spatially 
uniform steady state distribution only up to a critical total population density. (2) Beyond this critical 
density, the uniform distribution can not be maintained, and massively dense gregarious clusters form. 
(3) Linear stability analysis allows us to understand how the critical density depends on dimensionlcss 
ratios of the biological parameters. This dependence is summarized in Fig. [2] Our analysis also yields 
the most unstable cluster spacing (from the wave number of the most unstable modes). (4) Numerical 
simulations illustrate the rapid transitions that take place once gregarization is initiated. Dense packs 
of gregarious locusts form and grow, and these move and sweep up solitarious locusts in their vicinity. 
(5) Via bulk approximation, we find estimates for the long-time mass fraction dynamics of solitarious 
and gregarious locusts. In the large population limit, the entire population will become gregarious. Bulk 
theory and simulations agree well, as shown in Fig.[SJ (6) Our model displays population-level hysteresis, 
via which the critical density at which a gregarious aggregation forms from a dispersed population can 
be significantly higher than the density at which a gregarious aggregation would break up, as shown in 

Fig.m 

Our results shed light on locust control strategies in two ways. First, given the mass gregarization 
that takes place past the point of linear instability, the density threshold for this instability is a crucial 
quantity. In accordance with the idea proposed in [42], our work identifies a threshold below which 
populations should be kept in order to avoid a gregarious outbreak (assuming biological parameters 
are known to a sufficiently accurate degree). Furthermore, we have shown how this population-level 
property depends on individual-level parameters, finding a nontrivial relationship. Second, the apparent 
population-level hysteresis shows that dispersing a gregarized band, perhaps by killing individuals with 
pesticides, is harder than preventing group formation in the first place in that band annihilation requires 
a significantly lower locust density. In short, hysteresis implies that prevention could be more easily 
achieved than control. 

Like all models, ours has its limitations. We did not include features of the environment such as vege- 
tation, shown to have important influence on local crowding and hence gregarization. Our simplifications 
lead to mathematical tractability, while limiting the direct biological relevance of the model at present. In 
the field, locusts encounter patchy vegetation and other environmental influences, and adding such factors 
to the model would make it more relevant to field experiments. Since we have not explicitly included 
resource gradients or other environmental cues, we do not here recapitulate the long-range motion of 
locust bands, but merely their formation and clustering. Including environmental factors constitutes an 
extension of the current framework. Similarly, simulations in two spatial dimensions are more challenging 
and remain open for future investigation. 

Our work suggests several future biological experiments. First, as always, more accurate knowledge 
of model inputs would lead to better results. For our model, key inputs include the social interaction 
parameters, namely the length scales (r s , r g , and a g ) and interaction amplitudes (R s , R g , and A g ) in 
Eqn. ([5]) that we inferred from careful experiments such as those in [21] ■ However, to our knowledge, most 
of these parameters have not been directly measured in experiments on individuals. Second, we encourage 
observations of macroscopic group properties that could be compared to outputs of our model. These 
outputs include densities and sizes of bands. Additional quantitative field measurements along the lines 
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of [43] could help validate and refine our model. Finally, we can imagine experiments that would probe 
important aspects of the system dynamics (as opposed to physical properties of the bands themselves). 
Hopper bands arc known to undergo complicated dynamics, including splitting and merging [3J. BBC 
video shows an example of such phenomena in Locustana pardalina bands [44] . More accurate data 
for the dynamics of wild groups, including times for group formation and distances between merging 
bands and tributaries, could be compared to clumping time and length scales identified by our model. 
We are especially curious about experiments in which the critical average density for population-level 
gregarization and clumping might be probed in a controlled lab experiment, perhaps by slowly adding 
solitarious individuals into a large arena. Experimental measurements like those we have mentioned here 
would also motivate future two-dimensional extensions of our model where the streaming dynamics of 
hopper bands, the effects of the environment, and other stimuli could be more fully explored. 
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Figure 1. Spatially homogeneous steady states (HSS). (A) Spatially homogeneous solitarious 
(so, solid blue) and gregarious (go-, broken green) steady state locust density as the total locust density 
(po) is varied. For each set of curves, the middle curve represents the solution for our default phase 
change parameters, k\ — = 65 locusts/m? and Si = 62 = 0.25 hr~ l . The bottom and top curve in 
each set show parameter sensitivity; they are the 25th and 75th percentile values for so and go over 
10,000 parameter sets of {Si, 82, k\, k2} sampled from uniform distributions centered at the default 
values and varying by ±30%. In both the thin grey and red regions, the HSS is linearly unstable to small 
perturbations. Additionally, in the red region, go > sq, while in the grey and white regions, the opposite 
holds. (B) A blow-up of the boxed transition region in (A) around which the value of g overtakes s . 
The dashed black vertical lines indicate the 25th and 75th percentile for this transition. The solid 
purple vertical lines indicate the 25th and 75th percentile values for the onset of linear instability. 
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Figure 2. Linear stability of spatially homogeneous steady state (HSS) solutions. The 

dimensionless phase change parameter ratios 7 = 81/82 and K = fci/fca vary along the horizontal and 
vertical (we have used log axes). The contours indicate the critical value of rescaled density tp* = pl/k 2 - 
For rescaled densities greater than that value, the HSS solution is unstable. The critical rescaled 
density is monotonically increasing in both 7 and K . The arrow along the horizontal axis indicates the 
direction 7 moves if the relative rate of gregarization is increased (faster gregarization) . The arrow 
along the vertical axis indicates the direction K moves if the relative density threshold for gregarization 
is decreased (easier gregarization). For an accurate biological interpretation, one must multiply ip* by 
&2 in order to obtain the unsealed critical density p^. 



Figure 3. Maximally unstable perturbation wave number q max for homogeneous steady 
states with total density pq. Similar to Fig. [I] the middle, bottom, and top curves show results for 
the 25th and 75th percentile as computed from 10,000 random parameter draws centered around our 
default parameter set. At low densities, there are no unstable perturbation wave numbers. Just past 
the critical density p^, q nia , x increases rapidly and then plateaus. For our default parameters, q ma , x 
asymptotes to 8.89 mT 1 corresponding to a length scale of 0.71 m. 
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Figure 4. Numerical simulations of Eqs. Snapshots depict the numerical solution of 
Eqs. ©-([7]) at different times t (in hr) on a periodic domain of length L = 3 m with the default set of 
phase change parameters. See also Fig. [5] for a comparison with the alternative parameter set. The 
solitarious (gregarious) density (in locusts / m) as a function of spatial position (in m) is shown in blue 
(green). The total population mass is M = 50 locusts and the initial condition is set at g(x,t = 0) = 
and s(x, t = 0) given by a random perturbation centered around s = M/L. The top row of panels shows 
the fast smoothing of the initial state, and the subsequent evolution. Gregarization (approximately) 
occurs according to the spatially homogeneous version of Eqn. © , as can be seen up until the second 
row of panels, where the small instability becomes significant. Two compactly supported clumps of 
gregarious locusts form, superposed on a very sparse population of solitarious insects. In the third row, 
the gregarious group travels as a propagating pulse, and eventually stops. During this stage, the 
gregarious and solitarious populations are essentially non-overlapping in space. As shown in Fig.|6l the 
group continues to slowly gregarizc after it becomes stationary. 
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Figure 5. Numerical simulations of Eqs. ([3|). Similar to Fig. 01 snapshots at different times t (in 
hours), but for the alternative set of phase change parameters. Note that three, rather than two clumps 
of gregarious locusts form at intermediate times. This simulation is continued until t = 80 hr (last 
frame) to show the stability of the final cluster of gregarious locusts. 
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Figure 6. Population-level phase change over time. Mass fractions 4is,4>g of solitarious (blue) 
and gregarious (green) locusts as a function of time (in hours) for the numerical simulations of Fig.0] 
and Fig. [5] (A) Default set of phase change parameters, corresponding to the simulation in Fig.|4j (B) 
Alternative set of phase change parameters, corresponding to the simulation in Fig. [5] For both cases, 
at early times, these mass dynamics are well-approximated by the spatially homogeneous version of the 
governing equations Eqs. (j3|), whose solution, Eqn. (|T8|) . is shown as dotted curves. At late times, the 
mass dynamics are approximately described by the spatially segregated bulk theory of Eqn. (|19[) . whose 
solution is shown as dashed curves. 
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Figure 7. Population-level hysteresis as a function of average density po. Gregarious mass 
fraction <fi g as the average density po (total mass M divided by domain length L) is varied as a control 
parameter. We use our alternative set of phase change and social interaction parameters, as in Fig.[5j 
The solid (dotted) red curve represents the stable (unstable) homogeneous steady state solution, as 
calculated via linear stability analysis. As p a passes through the point of linear instability (marked with 
an asterisk) the solution jumps up to the green curve, which represents compactly supported gregarious 
aggregations obtained via numerical simulation, similar to the final states of Fig. U and Fig.[5j As p a is 
decreased by slowly subtracting mass from aggregations on the green curve, the system remains on the 
upper branch even for values of po sufficiently small as to be in the regime where the uniform state is 
stable, thus demonstrating dynamical population-level hysteresis. 



